In [2]:
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as mcolors
import pygrib as pb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import glob

#Get data path
datapath = glob.glob("./21072000/M-*grb2")
print(datapath)
['./21072000\\M-A0064-21072000-000.grb2', './21072000\\M-A0064-21072000-006.grb2', './21072000\\M-A0064-21072000-012.grb2', './21072000\\M-A0064-21072000-018.grb2', './21072000\\M-A0064-21072000-024.grb2', './21072000\\M-A0064-21072000-030.grb2', './21072000\\M-A0064-21072000-036.grb2', './21072000\\M-A0064-21072000-042.grb2', './21072000\\M-A0064-21072000-048.grb2', './21072000\\M-A0064-21072000-054.grb2', './21072000\\M-A0064-21072000-060.grb2', './21072000\\M-A0064-21072000-066.grb2', './21072000\\M-A0064-21072000-072.grb2', './21072000\\M-A0064-21072000-078.grb2', './21072000\\M-A0064-21072000-084.grb2']
In [10]:
"""
利用basemap套件繪製圖形,在read_grib範例中,可以得知變數的位置,這邊畫雨量跟10米風,初始場沒有雨量值
需要的步驟,
1.讀取grib資料,取出所需資料,包含2米溫度、10米u風、10米v風、經緯度
2.投影相關設定值(lon_0, lat_0, lat_1, lat_2)
  lon_0的key:LoVInDegrees
  lat_0的key:LaDInDegrees
  lat_1的key:Latin1InDegrees
  lat_2的key:Latin2InDegrees
3.繪製圖形
  這邊會先設定雨量值的colorbar
"""
cwb_data=['None','#9BFFFF','#00CFFF','#0198FF','#0165FF','#309901','#32FF00','#F8FF00','#FFCB00','#FF9A00',\
          '#FA0300','#CC0003','#A00000','#98009A','#C304CC','#F805F3','#FECBFF']
clevs = [0,1,2,6,10,15,20,30,40,50,70,90,110,130,150,200,300,400]
cmap = mcolors.ListedColormap(cwb_data,'preecipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)

cnt = 0 # 只是為了抓取1測的固定值用 
for grb in range(len(datapath)-1):
    if cnt == 0:
        forget = pb.open(datapath[grb]).select()[:][0]
        analDate = forget.analDate.strftime("%y%m%d%H")
        lon_0, lat_0 = forget["LoVInDegrees"], forget["LaDInDegrees"]
        lat_1, lat_2 = forget["Latin1InDegrees"], forget["Latin2InDegrees"]
        lats, lons = forget.latlons()
        proj = Basemap(projection="lcc",resolution='l',lat_0 = lat_0 ,lon_0 = lon_0 ,lat_1 = lat_1, \
              lat_2 = lat_2, llcrnrlat = lats[0,0], llcrnrlon = lons[0,0], urcrnrlat = lats[-1,-1], urcrnrlon = lons[-1,-1])
        cnt = 1
    basicvar1 = pb.open(datapath[grb]).select()[:]
    basicvar2 = pb.open(datapath[grb+1]).select()[:]
    fcst = str(basicvar2[0].forecastTime)
    valDate = basicvar2[0].validDate.strftime("%y%m%d%H")
    accrain = basicvar2[61]["values"] - basicvar1[61]["values"]
    u10m_data = basicvar2[66]["values"]
    v10m_data = basicvar2[67]["values"]
    fig = plt.figure(figsize=(16,12))
    proj.drawcoastlines(linewidth=0.9, color='slategrey')
    meridians = np.arange(60.,170.,5.)#設定繪製經度線的起點、終點及間隔
    parallels = np.arange(-10.,60.,5.)#設定繪製緯度線的起點、終點及間隔
    proj.drawmeridians(meridians, labels = [0,0,0,1],color='grey',linewidth=0.8)#繪製經度線 label = [左,右,上,下]
    proj.drawparallels(parallels, labels = [1,0,0,0],color='grey',linewidth=0.8)#繪製緯度線
    ctf = proj.contourf(lons, lats, accrain, clevs[1:len(clevs)-1], cmap=cmap, norm=norm, extend='both',latlon=True)
    gap = 20 #沒有用gap的話,會很慘
    proj.barbs(lons[::gap,::gap], lats[::gap,::gap], u10m_data[::gap,::gap], v10m_data[::gap,::gap] ,\
               length=5,barbcolor="navy", latlon=True)
    cbar = plt.colorbar(ctf,orientation='horizontal',pad=0.05,fraction=0.05)
    cbar.set_label("$^\circ$C",size=14)
    plt.title("Initial Time: {}".format(analDate), loc="left",fontsize=16)
    plt.title("Forecast Time:{}".format(fcst),loc="center",fontsize=16)
    plt.title("Valid Time:{}".format(valDate),loc="right",fontsize=16)
    plt.suptitle("In-Fa Typhoon",y=0.85,fontsize=20)
    plt.show()
    plt.clf()
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>